%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reproduces Appendix Figure A.2 Labor Tax Shock
%
% Cloyne, Dimsdale and Hürtgen (2024)
% "Are Tax Cuts Contractionary at the Zero Lower Bound?
%  Evidence from a Century of Data"
%
% James Cloyne, Nicholas Dimsdale and Patrick Hürtgen, 2024
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

%% set path and solve model using dynare and Occbin

set(0,'DefaultLineLineWidth',2)
addpath('c:/dynare/5.5/matlab')

dynare Occbin_NK_baseline_LaborTax console

%% generate draws from prior distribution

draw2=2000;         % number of actual draws
draws=2100;         % extra number of draws in case some draws are not valid

draw_i=1;
go_flag=0;
go = true;

[beta1,beta2]=betaAB(0.5,0.2); 
[beta3,beta4]=betaAB(0.6,0.15);
[gam1,gam2]  =gammaAB(2,0.5);
[beta5,beta6]=betaAB(0.8,0.1); 

while draw_i < draws && go

    display(['Draw ', num2str(draw_i)]);

    omegap=random('Beta',beta1,beta2);

    phipi=random('Normal',1.5,0.2); 
    while phipi<1.05
        phipi=random('Normal',1.5,0.2); 
    end

    phiy=random('Normal',0.125,0.05); 
    while phiy<0
        phiy=random('Normal',0.125,0.05); 
    end

    rhot1=random('Beta',beta5,beta6);

    gamma  = random('Normal',1.5,0.25);
    adj    = random('Normal',6,1.5);
    xi     = random('gamma',gam1,gam2);

    theta=0;
    chip=0;
    psi=1;

    actualdraws.omegap(draw_i)=omegap;
    actualdraws.phipi(draw_i)=phipi;
    actualdraws.phiy(draw_i)=phiy;
    actualdraws.gamma(draw_i)=gamma;
    actualdraws.xi(draw_i)=xi;
    actualdraws.adj(draw_i)=adj;
    actualdraws.rhot1(draw_i)=rhot1;

    draw_i=draw_i+1;
end

%% start priod predictive loop

tic
draw_i=1;
draws_skipped=[];
iter2=1;

while draw_i< draws

    drop_draw=0;

    omegap= actualdraws.omegap(draw_i);
    phiy  = actualdraws.phiy(draw_i);
    phipi = actualdraws.phipi(draw_i);
    xi    = actualdraws.xi(draw_i);
    gamma = actualdraws.gamma(draw_i);
    adj  = actualdraws.adj(draw_i);
    rhot1= actualdraws.rhot1(draw_i);

    set_param_value('rhot1' ,rhot1);
    set_param_value('omegap' ,omegap);
    set_param_value('phipi' ,phipi);
    set_param_value('phiy' ,phiy);
    set_param_value('adj' ,adj);
    set_param_value('xi' ,xi);
    set_param_value('gamma' ,gamma);

    options_.occbin.simul.SHOCKS=shock_mat.matrix_1;

    [oo_, out]  = occbin.solver(M_, oo_, options_);
    if out.error_flag %solution was not successful
        drop_draw=1;
    end

    options_.occbin.simul.SHOCKS=shock_mat.matrix_2;

    [oo3_, out]  = occbin.solver(M_, oo_, options_);
    if out.error_flag %solution was not successful
        drop_draw=1;
    end

    if drop_draw==0

        IRF.linear_Y=oo_.occbin.simul.linear(1:50+0,strmatch('y',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('y',M_.endo_names,'exact'));
        IRF.piecewise_Y=oo_.occbin.simul.piecewise(1:50+0,strmatch('y',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('y',M_.endo_names,'exact'));

        IRF.linear_INFC=oo_.occbin.simul.linear(1:50+0,strmatch('infc',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('infc',M_.endo_names,'exact'));
        IRF.piecewise_INFC=oo_.occbin.simul.piecewise(1:50+0,strmatch('infc',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('infc',M_.endo_names,'exact'));

        IRF.linear_RR=oo_.occbin.simul.linear(1:50+0,strmatch('rr',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('rr',M_.endo_names,'exact'));
        IRF.piecewise_RR=oo_.occbin.simul.piecewise(1:50+0,strmatch('rr',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('rr',M_.endo_names,'exact'));

        Simulation.y1(:,draw_i) = 100*IRF.linear_Y;
        Simulation.y2(:,draw_i) = 100*IRF.piecewise_Y;

        Simulation.rr1(:,draw_i) = 400*IRF.linear_RR;
        Simulation.rr2(:,draw_i) = 400*IRF.piecewise_RR;

        Simulation.invfc1(:,draw_i) = 100*IRF.linear_INFC;
        Simulation.invfc2(:,draw_i) = 100*IRF.piecewise_INFC;

    elseif drop_draw==1
        draws_skipped(iter2)=draw_i;
        iter2=iter2+1;
        drop_draw=0;
    end
    draw_i=draw_i + 1
end

toc

%% generate figures

start1=5;
end1=5+12;

xmin=0;
xmax=13;
stepsize=1;

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(1);
fanChart(0:size(Simulation.y2(start1:end1,:),1)-1, Simulation.y2(start1:end1,:)); hold on;
%title('GDP at ZLB')
set(gca,'TickDir','in'); 
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA2_GDP_ZLB_LaborTaxes');
print(filepath1, '-dpdf', '-fillpage'); 
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(2);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.invfc2(start1:end1,:));hold on;
%title('Inflation at the ZLB')
set(gca,'TickDir','in'); 
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA2_INFL_ZLB_LaborTaxes');
print(filepath1, '-dpdf', '-fillpage'); 
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(3);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.rr2(start1:end1,:));hold on;
%title('Real rate at ZLB')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
set(gca,'TickDir','in'); 
xlabel('Quarters')
set(gca,'YLIM', [-10 70]);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA2_RR_ZLB_LaborTaxes');
print(filepath1, '-dpdf', '-fillpage'); 
print(filepath1, '-depsc'); 







